# pacmanパッケージを利用しない場合
# library(tidyverse)
# library(MatchIt)
# ...
pacman::p_load(tidyverse,
MatchIt,
WeightIt,
BalanceR,
cobalt,
summarytools,
modelsummary,
kableExtra)
# cobaltパッケージが提供するデータセットの読み込み
data("lalonde", package = "cobalt")
lalonde_df <- lalondeマッチングにおける古典的なデータセット、lalondeを読み込みました。data(lalonde, package = "cobalt")を入力するだけで、{cobalt}パッケージ内のlaondeという名前のデータフレームが作業環境上に読み込まれます1。このデータをDay2_dfという名のオブジェクトとして改めて保存しておきます。
それでは、データの中身を確認してみましょう。普通にコンソール上でlalondeを打っても良いですが、DTパッケージのdatatable()を使うと見やすい表形式で出力されます。
# 今後、datatable()関数を使う予定はないので、DTパッケージを読み込まず、
# DT::datatable()を使用する
DT::datatable(lalonde_df)分析に入る前に、名目変数である人種(race)をダミー変数にします。raceは3種類の値で構成されているため、生成するダミー変数も3つとなります。
lalonde_df <- lalonde_df %>%
mutate(white = if_else(race == "white", 1, 0),
black = if_else(race == "black", 1, 0),
hispanic = if_else(race == "hispan", 1, 0),
# 作成された以上の3変数をraceの後へ
.after = race)続いて、{summarytools}記述統計量を確認してみましょう。Rコンソール上で確認するならdescr()を、R Markdown内に埋め込むならdfSummary()が良いでしょう。
# descr()を使う場合
descr(lalonde_df)## Descriptive Statistics
## lalonde_df
## N: 614
##
## age black educ hispanic married nodegree re74 re75
## ----------------- -------- -------- -------- ---------- --------- ---------- ---------- ----------
## Mean 27.36 0.40 10.27 0.12 0.42 0.63 4557.55 2184.94
## Std.Dev 9.88 0.49 2.63 0.32 0.49 0.48 6477.96 3295.68
## Min 16.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00
## Q1 20.00 0.00 9.00 0.00 0.00 0.00 0.00 0.00
## Median 25.00 0.00 11.00 0.00 0.00 1.00 1042.33 601.55
## Q3 32.00 1.00 12.00 0.00 1.00 1.00 7891.93 3254.81
## Max 55.00 1.00 18.00 1.00 1.00 1.00 35040.07 25142.24
## MAD 8.90 0.00 1.48 0.00 0.00 0.00 1545.36 891.86
## IQR 12.00 1.00 3.00 0.00 1.00 1.00 7888.50 3248.99
## CV 0.36 1.24 0.26 2.75 1.19 0.77 1.42 1.51
## Skewness 1.02 0.43 -0.72 2.37 0.34 -0.54 1.58 2.29
## SE.Skewness 0.10 0.10 0.10 0.10 0.10 0.10 0.10 0.10
## Kurtosis 0.14 -1.82 1.62 3.64 -1.89 -1.71 1.95 7.13
## N.Valid 614.00 614.00 614.00 614.00 614.00 614.00 614.00 614.00
## Pct.Valid 100.00 100.00 100.00 100.00 100.00 100.00 100.00 100.00
##
## Table: Table continues below
##
##
##
## re78 treat white
## ----------------- ---------- -------- --------
## Mean 6792.83 0.30 0.49
## Std.Dev 7470.73 0.46 0.50
## Min 0.00 0.00 0.00
## Q1 237.91 0.00 0.00
## Median 4759.02 0.00 0.00
## Q3 10923.35 1.00 1.00
## Max 60307.93 1.00 1.00
## MAD 7055.72 0.00 0.00
## IQR 10655.31 1.00 1.00
## CV 1.10 1.52 1.03
## Skewness 1.55 0.86 0.05
## SE.Skewness 0.10 0.10 0.10
## Kurtosis 4.33 -1.26 -2.00
## N.Valid 614.00 614.00 614.00
## Pct.Valid 100.00 100.00 100.00
# dfSummary()はHTML形式で出力されるため、R Markdownに埋め込む
dfSummary(lalonde_df, style = "grid",
plain.ascii = FALSE, graph.magnif = 0.85) %>%
print(method = "render", headings = FALSE) | No | Variable | Stats / Values | Freqs (% of Valid) | Graph | Valid | Missing | ||||||||||||
|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|---|
| 1 | treat [integer] | Min : 0 Mean : 0.3 Max : 1 |
|
614 (100.0%) | 0 (0.0%) | |||||||||||||
| 2 | age [integer] | Mean (sd) : 27.4 (9.9) min < med < max: 16 < 25 < 55 IQR (CV) : 12 (0.4) | 40 distinct values | 614 (100.0%) | 0 (0.0%) | |||||||||||||
| 3 | educ [integer] | Mean (sd) : 10.3 (2.6) min < med < max: 0 < 11 < 18 IQR (CV) : 3 (0.3) | 19 distinct values | 614 (100.0%) | 0 (0.0%) | |||||||||||||
| 4 | race [factor] | 1. black 2. hispan 3. white |
|
614 (100.0%) | 0 (0.0%) | |||||||||||||
| 5 | white [numeric] | Min : 0 Mean : 0.5 Max : 1 |
|
614 (100.0%) | 0 (0.0%) | |||||||||||||
| 6 | black [numeric] | Min : 0 Mean : 0.4 Max : 1 |
|
614 (100.0%) | 0 (0.0%) | |||||||||||||
| 7 | hispanic [numeric] | Min : 0 Mean : 0.1 Max : 1 |
|
614 (100.0%) | 0 (0.0%) | |||||||||||||
| 8 | married [integer] | Min : 0 Mean : 0.4 Max : 1 |
|
614 (100.0%) | 0 (0.0%) | |||||||||||||
| 9 | nodegree [integer] | Min : 0 Mean : 0.6 Max : 1 |
|
614 (100.0%) | 0 (0.0%) | |||||||||||||
| 10 | re74 [numeric] | Mean (sd) : 4557.5 (6478) min < med < max: 0 < 1042.3 < 35040.1 IQR (CV) : 7888.5 (1.4) | 358 distinct values | 614 (100.0%) | 0 (0.0%) | |||||||||||||
| 11 | re75 [numeric] | Mean (sd) : 2184.9 (3295.7) min < med < max: 0 < 601.5 < 25142.2 IQR (CV) : 3249 (1.5) | 356 distinct values | 614 (100.0%) | 0 (0.0%) | |||||||||||||
| 12 | re78 [numeric] | Mean (sd) : 6792.8 (7470.7) min < med < max: 0 < 4759 < 60307.9 IQR (CV) : 10655.3 (1.1) | 457 distinct values | 614 (100.0%) | 0 (0.0%) |
Generated by summarytools 0.9.9 (R version 4.1.0)
2021-07-27
treat変数の値ごとに結果変数であるre78の平均値を計算してみましょう。
Diff_Mean_df <- lalonde_df %>%
group_by(treat) %>%
summarise(Outcome = mean(re78),
.groups = "drop")
Diff_Mean_df## # A tibble: 2 x 2
## treat Outcome
## <int> <dbl>
## 1 0 6984.
## 2 1 6349.
この結果を可視化するコードは以下の通りです。
Diff_Mean_df %>%
ggplot() +
geom_bar(aes(x = treat, y = Outcome),
stat = "identity", width = 0.5) +
geom_label(aes(x = treat, y = Outcome,
label = round(Outcome, 3))) +
labs(x = "Treatment",
y = "Outcome (US Dollars)") +
# scale_x_continuous()を使って0/1をControl/Treatmentに置換する
# 目盛りはX軸上の0と1、各目盛りのラベルはControlとTreatmentに
scale_x_continuous(breaks = c(0, 1), labels = c("Control", "Treatment")) +
coord_cartesian(xlim = c(-0.5, 1.5))treat == 0の回答者、つまり職業訓練を受けていない回答者の平均所得は約6984ドル、treat == 1の回答者、つまり職業訓練を受けた回答者の平均所得は約6394ドルです。その差は約-650ドルですが、職業訓練を受けた回答者の方が低所得になっています。これは直感的に納得できる結果ではないでしょう。むろん、実際、職業訓練が所得を減らす可能性もありますが、今回の結果はより詳しく分析してみる価値はあるでしょう。
ちなみに、以上の結果は単回帰分析からも確認可能です (ただし、統計的に有意ではない)。
Raw_Fit1 <- lm(re78 ~ treat, data = lalonde_df)
# modelsummaryを使った推定結果の要約
# デフォルトでは点推定値と標準誤差が出力されるが、
# ここでは点推定値と95%信頼区間を出力する。
modelsummary(Raw_Fit1,
statistic = "[{conf.low}, {conf.high}]",
conf_level = 0.95)| Model 1 | |
|---|---|
| (Intercept) | 6984.170 |
| [6275.791, 7692.549] | |
| treat | -635.026 |
| [-1925.544, 655.492] | |
| Num.Obs. | 614 |
| R2 | 0.002 |
| R2 Adj. | 0.000 |
| AIC | 12698.7 |
| BIC | 12712.0 |
| Log.Lik. | -6346.371 |
| F | 0.934 |
この直感的でない結果は、もしかしたらセレクションバイアスが原因かも知れません。職業訓練の対象が元々非常に所得が低い被験者になっている可能性があります。たとえば、下の図のように職業訓練の有無が教育水準や人種、これまでの所得などと関係しているとしましょう。これらの要因は回答者の現在所得にも関係していると考えられます。この場合、処置有無と所得の間には内生性が存在することになります。
本当にそうなのか、共変量のバランスチェックをしてみましょう。もし、処置有無によって回答者の社会経済的要因に大きな差があれば、内生性が存在する証拠になります。ここでは誰かが作成しました{BalanceR}パッケージを使います。
Blc_Chk <- lalonde_df %>%
BalanceR(group = treat, cov = c(age, educ, white:re75)){BalanceR}パッケージで共変量を指定する際、:演算子が使用可能です。white:re75は、データセットのwhiteからre75変数までをすべて指定することを意味します。Rコンソール上で変数名のみを出力するにはnames()関数を使用します。
names(lalonde_df)## [1] "treat" "age" "educ" "race" "white" "black"
## [7] "hispanic" "married" "nodegree" "re74" "re75" "re78"
whiteからre75までblack、hispanic、marriedなどの変数がありますが、これらを一々指定せず、:を使えば一気に指定することができます。それではバランスチェックの結果を確認してみましょう。
Blc_Chk## Covariate Mean:0 SD:0 Mean:1 SD:1 SB:0-1
## 1 age 28.030 10.787 25.816 7.155 24.190
## 2 educ 10.235 2.855 10.346 2.011 -4.476
## 3 white 0.655 0.476 0.097 0.297 140.799
## 4 black 0.203 0.403 0.843 0.365 -167.083
## 5 hispanic 0.142 0.350 0.059 0.237 27.740
## 6 married 0.513 0.500 0.189 0.393 72.076
## 7 nodegree 0.597 0.491 0.708 0.456 -23.549
## 8 re74 5619.237 6788.751 2095.574 4886.620 59.575
## 9 re75 2466.484 3291.996 1532.055 3219.251 28.700
| Covariate | Mean:0 | SD:0 | Mean:1 | SD:1 | SB:0-1 |
|---|---|---|---|---|---|
| age | 28.030 | 10.787 | 25.816 | 7.155 | 24.190 |
| educ | 10.235 | 2.855 | 10.346 | 2.011 | -4.476 |
| white | 0.655 | 0.476 | 0.097 | 0.297 | 140.799 |
| black | 0.203 | 0.403 | 0.843 | 0.365 | -167.083 |
| hispanic | 0.142 | 0.350 | 0.059 | 0.237 | 27.740 |
| married | 0.513 | 0.500 | 0.189 | 0.393 | 72.076 |
| nodegree | 0.597 | 0.491 | 0.708 | 0.456 | -23.549 |
| re74 | 5619.237 | 6788.751 | 2095.574 | 4886.620 | 59.575 |
| re75 | 2466.484 | 3291.996 | 1532.055 | 3219.251 | 28.700 |
いくつか怪しい箇所があります。たとえば、treat == 0の回答者において黒人の割合は約20%ですが、treat == 1のそれは約85%です。つまり、黒人ほどより職業訓練を受ける傾向があることを意味します。実際、標準化バイアスは-167という、非常に大きい数値を示しています。この結果を表としてまとめてみましょう。
# 標準化バイアスの可視化
# 絶対値変換。SB = 25に破線
plot(Blc_Chk, abs = TRUE, vline = 25)かなり緩めの基準である25を採用しても、人種(黒人)、結婚有無、74年の所得のバランスが非常によくありません。内生性があると判断して良いでしょう。以下では様々な方法を使い、この内生性に対処していきたいと思います。
つづいて、重回帰分析もしてみましょう。用いる共変量は年齢、教育水準、黒人ダミー、ヒスパニックダミー、既婚ダミー、学位なしダミー、74・75年の所得です。lm()関数を使用し、78年の所得をこちらの変数に回帰させます。
\[ \begin{align} \textrm{re78} = & \beta_0 + \beta_1 \textrm{treat} + \beta_2 \textrm{age} + \beta_3 \textrm{educ} + \\ & \beta_4 \textrm{black} + \beta_5 \textrm{hispanic} + \beta_6 \textrm{married} + \beta_7 \textrm{nodegree} + \beta_8 \textrm{re74} + \beta_9 \textrm{re75}. \end{align} \]
Raw_Fit2 <- lm(re78 ~ treat + age + educ + black + hispanic + married +
nodegree + re74 + re75, data = lalonde_df)
modelsummary(Raw_Fit2)| Model 1 | |
|---|---|
| (Intercept) | 66.515 |
| (2436.746) | |
| treat | 1548.244 |
| (781.279) | |
| age | 12.978 |
| (32.489) | |
| educ | 403.941 |
| (158.906) | |
| black | -1240.644 |
| (768.764) | |
| hispanic | 498.897 |
| (941.943) | |
| married | 406.621 |
| (695.472) | |
| nodegree | 259.817 |
| (847.442) | |
| re74 | 0.296 |
| (0.058) | |
| re75 | 0.232 |
| (0.105) | |
| Num.Obs. | 614 |
| R2 | 0.148 |
| R2 Adj. | 0.135 |
| AIC | 12617.5 |
| BIC | 12666.1 |
| Log.Lik. | -6297.752 |
| F | 11.636 |
共変量を統制したら処置変数の係数は約1548.244ドルです。先ほどの単回帰分析の結果とは違って、統計的に有意な正の効果が確認されていますね。
重回帰分析は非常にシンプルで便利な分析方法ですが、いくつかの欠点があります。まず、重回帰分析は変数間の関係(線形結合)および誤差項の分布(平均0の正規分布)などを仮定したパラメトリック分析になります。この場合、同じ共変量を持たないケースであっても、勝手に予測を行います。重回帰分析における処置変数の解釈は「他の共変量がすべて同じ」場合の処置効果になります。これは、共変量がすべて同じ場合における(最初に見た)単純差分のようなものです。しかし、「他の共変量がすべて同じ」ケースが存在しない可能性があります。特に、共変量が多く、連続変数の場合、共変量がすべて同じことは実質あり得ないか、非常に少ないケースに限定されます。一方、マッチングを行うと、「他の共変量がすべて同じ」、または「非常に似ている」ケース間で比較を行います。
ここでは以下の4つの手法を紹介します。
まずは、マハラノビス最近傍マッチングから始めましょう。だいたいのマッチング手法は{MatchIt}パッケージで解決できます。マッチングデータセットを作成する関数はmatchit()関数であり、使い方は
matchit(処置変数 ~ 共変量1 + ... + 共変量k,
data = データフレーム名, estimand = "ATT",
method = "nearest", distance = "mahalanobis")のように入力します。method = "nearest"は最近傍マッチングを意味し、distance = "mahalanobis"はマハラノビス距離を意味します。estimand = "ATT"はATTを推定することを意味します。{MatchIt}の最近傍マッチングの場合、"ATT"、または"ATC"のみ指定可能です。早速やってみましょう。
MH_Matching <- matchit(treat ~ age + educ + black + hispanic + married +
nodegree + re74 + re75, data = lalonde_df,
method = "nearest", distance = "mahalanobis",
estimand = "ATT")バランスチェックは数値的に確認する方法と視覚的にする方法があります。まずは、数値的に確認してみましょう。
# バランスチェック (Numerical)
summary(MH_Matching)##
## Call:
## matchit(formula = treat ~ age + educ + black + hispanic + married +
## nodegree + re74 + re75, data = lalonde_df, method = "nearest",
## distance = "mahalanobis", estimand = "ATT")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age 25.8162 28.0303 -0.3094 0.4400 0.0813
## educ 10.3459 10.2354 0.0550 0.4959 0.0347
## black 0.8432 0.2028 1.7615 . 0.6404
## hispanic 0.0595 0.1422 -0.3498 . 0.0827
## married 0.1892 0.5128 -0.8263 . 0.3236
## nodegree 0.7081 0.5967 0.2450 . 0.1114
## re74 2095.5737 5619.2365 -0.7211 0.5181 0.2248
## re75 1532.0553 2466.4844 -0.2903 0.9563 0.1342
## eCDF Max
## age 0.1577
## educ 0.1114
## black 0.6404
## hispanic 0.0827
## married 0.3236
## nodegree 0.1114
## re74 0.4470
## re75 0.2876
##
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age 25.8162 25.5351 0.0393 0.7451 0.0341
## educ 10.3459 10.4703 -0.0618 0.8918 0.0122
## black 0.8432 0.2595 1.6057 . 0.5838
## hispanic 0.0595 0.1297 -0.2971 . 0.0703
## married 0.1892 0.4486 -0.6625 . 0.2595
## nodegree 0.7081 0.6108 0.2140 . 0.0973
## re74 2095.5737 2545.7229 -0.0921 1.2602 0.0748
## re75 1532.0553 1645.4710 -0.0352 1.3141 0.0386
## eCDF Max Std. Pair Dist.
## age 0.1243 0.3188
## educ 0.0973 0.2285
## black 0.5838 1.6949
## hispanic 0.0703 0.6172
## married 0.2595 1.0765
## nodegree 0.0973 0.2378
## re74 0.3405 0.2010
## re75 0.1838 0.1586
##
## Percent Balance Improvement:
## Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
## age 87.3 64.2 58.1 21.2
## educ -12.5 83.7 64.8 12.6
## black 8.8 . 8.8 8.8
## hispanic 15.1 . 15.1 15.1
## married 19.8 . 19.8 19.8
## nodegree 12.6 . 12.6 12.6
## re74 87.2 64.8 66.7 23.8
## re75 87.9 -511.2 71.2 36.1
##
## Sample Sizes:
## Control Treated
## All 429 185
## Matched 185 185
## Unmatched 244 0
## Discarded 0 0
表が3つでますが、上から(1)マッチング前の各共変量の情報、(2)マッチング後の各共変量の情報、(3)バランス改善の指標となります。(3)から見る場合、hispanic変数は改善率が100%となっています。実際、マッチング前は統制群と処置群におけるhispanic変数の平均値の差分は-0.0827でしたが、マッチング後は0となり、完全にマッチングされたことが分かります。他にも多くの共変量においてバランスの改善が見られます。eudcはむしろ悪化しましたが、他の変数はかなり改善されているので、このままいきましょう。
つづいて視覚的な方法です。
# バランスチェック (Graphical): QQプロットによるバランスチェック
plot(MH_Matching)各ポイントが45度直線上に並ぶのが理想的なバランシングです。バランスが改善されていることが確認できます。
もっと直感的なバランスチェックの方法として{BalanceR}パッケージと同様、標準化差分を使う方法です。マッチング前後のバランスを同時に確認するには{cobalt}パッケージのlove.plotが便利です。
love.plot(MH_Matching, thresholds = 0.25, abs = TRUE)thresholds引数は垂直線(破線)の位置、absは標準化差分を絶対値で示すことを意味します。マッチング後の標準化差分(Adjusted; 赤い点)が0.25より左側に位置している場合、バランスしていると判断できます2。他にもマッチング後の標準化差分がマッチング前(Unadjusted; 青い点)より改善されるいるか否かも判断できます。今回の例だと、大幅にバランスが改善されました。blackはまだバランスが取れておりませんが、それでも大幅に改善されていることが分かります。
それではATTを推定してみましょう。推定方法としてはノンパラメトリックな方法とパラメトリック方法がありますが、結果は変わりません。ノンパラメトリックな方法はペアごとの差分を計算し、その平均値を求める方法ですが、マッチング済みのデータに対し、処置変数を結果変数を回帰させることも、結果的には同じことを行うことになります。したがって、もっと簡単なパラメトリック方法でATTを推定します。ノンパラメトリック方法によるATTの推定は本資料の付録を参照してください。
回帰分析を行うためにはデータが必要ですね。つまり、マッチングされないケースをデータから除去する必要があります。したがって、マッチングされた後のデータセットを抽出する必要があります。使う関数はmatch.data()関数です。
# マッチング後データを抽出
MH_Data <- match.data(MH_Matching)マッチングデータが取れたら、その中身を確認してみましょう。
dim(MH_Data)## [1] 370 14
MH_Dataデータのサイズは370行14列です。この370行には意味があります。それは処置群の大きさの2倍という点です。多くの場合、マッチングから計算される処置効果はATEではなく、ATTです。したがって、処置群のデータを100%活用し、統制群から共変量が最も近いケースを抽出&マッチングすることになります。だから、マッチング後のサンプルサイズは処置群のサイズの2倍になります。
それでは職業訓練のATTを推定します。方法は簡単です。このマッチング後のデータを用い、単回帰分析を行うのみです。
# 処置効果の確認 (Model-Based)
MH_Fit <- lm(re78 ~ treat, data = MH_Data)
modelsummary(MH_Fit)| Model 1 | |
|---|---|
| (Intercept) | 5842.100 |
| (528.119) | |
| treat | 507.043 |
| (746.873) | |
| Num.Obs. | 370 |
| R2 | 0.001 |
| R2 Adj. | -0.001 |
| AIC | 7624.8 |
| BIC | 7636.6 |
| Log.Lik. | -3809.419 |
| F | 0.461 |
処置効果は約507.043ドルであり、付録にあるノンパラメトリック推定値と一致していることが分かります。今回の結果は重回帰分析よりも推定値が低めであり、統計的に有意に職業訓練の効果があったとは言えないという結果が得られましたね。
ちなみに、{MatchIt}パッケージを使った最近傍マッチングのの結果は行う度に変化することがあります。なぜなら、{MatchIt}パッケージを使った最近傍マッチングの場合、処置群 (統制群)から一つのケースを選択し、最も近い統制群 (処置群)とマッチングさせます。マッチングされたケースは次のステップからはマッチング対象から除外されます3。また、1:1マッチングの場合4、同距離に複数のマッチング対象があると、ランダムに1つのみを選択します。最近傍マッチングを用いる際は、複数推定を行い、推定が安定するかを確認し、不安定な場合は他の手法を使うか、k-最近傍マッチングなどを使ってみましょう。
CEMはマハラノビス最近傍マッチング同様、matchit()関数を使います。ただし、cemパッケージをインストールしておく必要があります。まず、Exact Matching類の手法は距離を図る必要がないので、distance引数は不要です。method = ...引数を"nearest" (最近傍)から"cem"に替えるだけです。推定可能な処置効果はATE、ATT、ATCですが、ここではATTを推定してみましょう。
マッチングをしたらmatch.data()でマッチングされたデータを抽出します。
CEM_Matching <- matchit(treat ~ age + educ + black + hispanic + married +
nodegree + re74 + re75, data = lalonde_df,
method = "cem", estimand = "ATT")
summary(CEM_Matching)##
## Call:
## matchit(formula = treat ~ age + educ + black + hispanic + married +
## nodegree + re74 + re75, data = lalonde_df, method = "cem",
## estimand = "ATT")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age 25.8162 28.0303 -0.3094 0.4400 0.0813
## educ 10.3459 10.2354 0.0550 0.4959 0.0347
## black 0.8432 0.2028 1.7615 . 0.6404
## hispanic 0.0595 0.1422 -0.3498 . 0.0827
## married 0.1892 0.5128 -0.8263 . 0.3236
## nodegree 0.7081 0.5967 0.2450 . 0.1114
## re74 2095.5737 5619.2365 -0.7211 0.5181 0.2248
## re75 1532.0553 2466.4844 -0.2903 0.9563 0.1342
## eCDF Max
## age 0.1577
## educ 0.1114
## black 0.6404
## hispanic 0.0827
## married 0.3236
## nodegree 0.1114
## re74 0.4470
## re75 0.2876
##
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## age 21.6462 21.2936 0.0493 0.8909 0.0130
## educ 10.3692 10.2795 0.0446 0.8894 0.0100
## black 0.8615 0.8615 0.0000 . 0.0000
## hispanic 0.0308 0.0308 0.0000 . 0.0000
## married 0.0308 0.0308 0.0000 . 0.0000
## nodegree 0.6769 0.6769 0.0000 . 0.0000
## re74 489.0612 697.6115 -0.0427 1.3916 0.0435
## re75 362.4365 520.8378 -0.0492 0.8501 0.0403
## eCDF Max Std. Pair Dist.
## age 0.0897 0.1373
## educ 0.1000 0.1949
## black 0.0000 0.0000
## hispanic 0.0000 0.0000
## married 0.0000 0.0000
## nodegree 0.0000 0.0000
## re74 0.3118 0.0968
## re75 0.1598 0.1635
##
## Percent Balance Improvement:
## Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
## age 84.1 85.9 84.0 43.1
## educ 18.8 83.3 71.2 10.2
## black 100.0 . 100.0 100.0
## hispanic 100.0 . 100.0 100.0
## married 100.0 . 100.0 100.0
## nodegree 100.0 . 100.0 100.0
## re74 94.1 49.7 80.6 30.2
## re75 83.0 -263.5 70.0 44.4
##
## Sample Sizes:
## Control Treated
## All 429. 185
## Matched (ESS) 41.29 65
## Matched 75. 65
## Unmatched 354. 120
## Discarded 0. 0
CEMの場合、マッチングされないブロックは捨てられるため、マハラノビス距離最近傍マッチングよりもサンプルサイズが小さくなります。また、処置群と統制群のサンプルサイズも不均衡になります。マッチング結果を見ると、処置群からは65ケース、統制群からは75サンプルのみ残ることになりました。
つづいて、cobalt::love.plot()を使用して、バランスチェックを行います。
love.plot(CEM_Matching, thresholds = 0.25, abs = TRUE)CEMの場合、マハラノビス距離最近傍マッチングよりもバランスが大きく改善されたことが分かります。その理由は簡単です。最近傍マッチングの場合、最も近いケースであれば、どれほど離れていてもマッチングされます。一方、CEMは正確マッチングの一種であるため、ある程度離れているケースを捨ててしまうからです。結局は共変量が非常に近いケースのみを残すことになります。
それでは、match.data()関数を使ってマッチング後のデータを抽出してみましょう。
CEM_Data <- match.data(CEM_Matching)処置効果の推定の前に、抽出されたデータの中身を見てましょう。元のデータにはweightsという列はありませんでしたが、いきなりできました。実はこれが、CEMに使う重み変数です。自動的に計算してデータセットに追加してくれたわけです5。
CEM_Data処置効果の推定は重み付き回帰分析 (WLS)で行います。既存のlm()関数にweight = ...引数を追加し、データ内の重み列をしてするだけです。
# 処置効果の確認
# weight = Wが追加されていることを確認
CEM_Fit <- lm(re78 ~ treat, data = CEM_Data, weights = weights)
# 処置効果の確認
modelsummary(CEM_Fit)| Model 1 | |
|---|---|
| (Intercept) | 5265.785 |
| (850.457) | |
| treat | 1070.907 |
| (1248.129) | |
| Num.Obs. | 140 |
| R2 | 0.005 |
| R2 Adj. | -0.002 |
| AIC | 2924.7 |
| BIC | 2933.5 |
| Log.Lik. | -1459.329 |
| F | 0.736 |
今回、処置効果の推定値は約1070.907ドルです。
傾向スコアマッチングも、これまでのコマンドとほぼ同じです。マハラノビス最近傍マッチングのコマンドからdistance = ...引数を抜けば、傾向スコアマッチングができます6。また、replace引数にTRUEを指定します。これは一度マッチングしたケースを捨てずに、次のマッチングにも使うことを意味します。デフォルトはFALSEとなっており、この場合、一度マッチングされたケースは二度とマッチングされません。したがって、replace = TRUEを指定した方が共変量バランスがより改善されます。ただし、有効サンプルサイズ(Effective Sample Size; ESS)が小さくなり、精度が悪くなる点、場合によっては特殊な標準誤差7を使う必要があるといった欠点があります。これは最近傍マッチングに共通する問題ですが、ここではこれまでの結果と比較するために既定値のままにしましょう。
PS_Matching <- matchit(treat ~ age + educ + black + hispanic + married +
nodegree + re74 + re75, data = lalonde_df,
method = "nearest", estimand = "ATT")
summary(PS_Matching)##
## Call:
## matchit(formula = treat ~ age + educ + black + hispanic + married +
## nodegree + re74 + re75, data = lalonde_df, method = "nearest",
## estimand = "ATT")
##
## Summary of Balance for All Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5774 0.1822 1.7941 0.9211 0.3774
## age 25.8162 28.0303 -0.3094 0.4400 0.0813
## educ 10.3459 10.2354 0.0550 0.4959 0.0347
## black 0.8432 0.2028 1.7615 . 0.6404
## hispanic 0.0595 0.1422 -0.3498 . 0.0827
## married 0.1892 0.5128 -0.8263 . 0.3236
## nodegree 0.7081 0.5967 0.2450 . 0.1114
## re74 2095.5737 5619.2365 -0.7211 0.5181 0.2248
## re75 1532.0553 2466.4844 -0.2903 0.9563 0.1342
## eCDF Max
## distance 0.6444
## age 0.1577
## educ 0.1114
## black 0.6404
## hispanic 0.0827
## married 0.3236
## nodegree 0.1114
## re74 0.4470
## re75 0.2876
##
##
## Summary of Balance for Matched Data:
## Means Treated Means Control Std. Mean Diff. Var. Ratio eCDF Mean
## distance 0.5774 0.3629 0.9739 0.7566 0.1321
## age 25.8162 25.3027 0.0718 0.4568 0.0847
## educ 10.3459 10.6054 -0.1290 0.5721 0.0239
## black 0.8432 0.4703 1.0259 . 0.3730
## hispanic 0.0595 0.2162 -0.6629 . 0.1568
## married 0.1892 0.2108 -0.0552 . 0.0216
## nodegree 0.7081 0.6378 0.1546 . 0.0703
## re74 2095.5737 2342.1076 -0.0505 1.3289 0.0469
## re75 1532.0553 1614.7451 -0.0257 1.4956 0.0452
## eCDF Max Std. Pair Dist.
## distance 0.4216 0.9740
## age 0.2541 1.3938
## educ 0.0757 1.2474
## black 0.3730 1.0259
## hispanic 0.1568 1.0743
## married 0.0216 0.8281
## nodegree 0.0703 1.0106
## re74 0.2757 0.7965
## re75 0.2054 0.7381
##
## Percent Balance Improvement:
## Std. Mean Diff. Var. Ratio eCDF Mean eCDF Max
## distance 45.7 -239.6 65.0 34.6
## age 76.8 4.6 -4.2 -61.1
## educ -134.8 20.4 31.2 32.1
## black 41.8 . 41.8 41.8
## hispanic -89.5 . -89.5 -89.5
## married 93.3 . 93.3 93.3
## nodegree 36.9 . 36.9 36.9
## re74 93.0 56.8 79.1 38.3
## re75 91.2 -800.7 66.3 28.6
##
## Sample Sizes:
## Control Treated
## All 429 185
## Matched 185 185
## Unmatched 244 0
## Discarded 0 0
処置効果の推定の前に、バランスが取れているかどうかを確認してみましょう。これまではQQプロットのみでバランスチェックをしましたが、今回はヒストグラムとJitterプロットで確認します。
# 1. QQプロットによるバランスチェック
plot(PS_Matching)QQプロットを見ると、バランスが多少改善されていることが分かります。
# 2. ヒストグラムによるバランスチェック
plot(PS_Matching, type = "hist")ヒストグラムからもバランスが改善されたことが確認できますね。出来れば、ヒストグラムを重ねて比較したいんですが、残念です。あとでオーバーラップされたヒストグラムの作成方法について解説します。
# 3. 散布図によるバランスチェック
par(mfrow = c(1, 1))
plot(PS_Matching, type = "jitter")## [1] "To identify the units, use first mouse button; to stop, use second."
## integer(0)
散布図をみると、傾向スコア (\(e\))が小さいケースがマッチングから除外されたことが分かります。実際、処置群には傾向スコアが小さいケースが非常に少ないため、これは妥当な結果でしょう。
むろん、おなじみのcobalt::love.plot()でもバランスチェックができます。
love.plot(PS_Matching, thresholds = 0.25, abs = TRUE)それでは、ATT推定のためにマッチング後のデータを抽出します。
# 傾向スコアマッチング後のデータセットを抽出
PS_Data <- match.data(PS_Matching)傾向スコアを用いたATTをの推定もこれまでと同様、回帰分析を用います。普通の回帰分析ですが、用いるデータは傾向スコアでマッチングされたデータのみになります8。
# 処置効果の推定
PS_Fit <- lm(re78 ~ treat, data = PS_Data)
modelsummary(PS_Fit)| Model 1 | |
|---|---|
| (Intercept) | 5454.776 |
| (516.396) | |
| treat | 894.367 |
| (730.295) | |
| Num.Obs. | 370 |
| R2 | 0.004 |
| R2 Adj. | 0.001 |
| AIC | 7608.2 |
| BIC | 7620.0 |
| Log.Lik. | -3801.114 |
| F | 1.500 |
処置効果は約894.367ドルでした。
ちなみに、講義スライド47〜48ページ目の例を{MatchIt}パッケージを使うと以下のようなコードになります。講義スライドでは一度マッチングされたケースを再利用したため、replace = TRUEを指定し、ATT/ATCを推定する際にweights =引数を指定する必要があります。
# 講義スライドのデータ
PS_df2 <- read.csv("Data/Day2_Data2.csv")
p47_mat <- matchit(Treat ~ Age + Edu, estimand = "ATT",
data = PS_df2, method = "nearest", replace = TRUE)
p48_mat <- matchit(Treat ~ Age + Edu, estimand = "ATC",
data = PS_df2, method = "nearest", replace = TRUE)
p47_data <- match.data(p47_mat)
p48_data <- match.data(p48_mat)
p47_fit <- lm(Outcome ~ Treat, data = p47_data, weights = weights)
p48_fit <- lm(Outcome ~ Treat, data = p48_data, weights = weights)
modelsummary(list("ATT" = p47_fit, "ATC" = p48_fit))| ATT | ATC | |
|---|---|---|
| (Intercept) | 4.909 | 4.308 |
| (0.830) | (0.599) | |
| Treat | 1.818 | 3.923 |
| (1.031) | (1.012) | |
| Num.Obs. | 17 | 20 |
| R2 | 0.172 | 0.455 |
| R2 Adj. | 0.116 | 0.425 |
| AIC | 76.7 | 93.8 |
| BIC | 79.2 | 96.8 |
| Log.Lik. | -35.343 | -43.906 |
| F | 3.108 | 15.031 |
手計算から得られた結果と一致することが分かります。
最後にATTでなく、ATEが推定の対象となるIPW推定量を計算してみましょう。{WeightIt}パッケージを使用すれば便利ですが、ここではあえて使わずにIPW推定量を計算してみましょう。こうすることでIPWが具体的にどのような仕組みで推定されるかが理解できるはずです。一通りの推定が終わりましたら{WeightIt}パッケージの使い方も紹介します。
まずは、傾向スコアの算出ですが、ここではロジスティック回帰分析で傾向スコアを計算します9。
# Logitで傾向スコアを推定
PS <- glm(treat ~ age + educ + black + hispanic + married + nodegree +
re74 + re75, data = lalonde_df,
family = binomial("logit"))続いて、既存のデータフレームに傾向スコアが格納された列を追加します。まず、既存のlalondeデータセットをMatching.dfという名前で複製します。
新しいデータフレームにPSという列を作成し、ここに傾向スコアを格納します。傾向スコアはロジスティック回帰分析のオブジェクトに$fitted.valuesを加えるだけで抽出できます。
# lalondeデータフレームを複製
ipw_df <- lalonde_df
# 傾向スコア推定値をデータフレーム内に格納
ipw_df$PS <- PS$fitted.values次は、各ケースの重みを計算し、Wという新しい列として追加します。傾向スコアは、もしケースが処置群 (treat == 1)なら、\(\frac{1}{e}\)を、統制群 (treat == 0)なら\(\frac{1}{(1-e)}\)となります。
# 重み変数を作成
## ケースが処置群なら、重み = 傾向スコアの逆数
## ケースが統制群あら、重み = (1 - 傾向スコア)の逆数
ipw_df <- ipw_df %>%
mutate(W = ifelse(treat == 1, (1 / PS), (1 / (1 - PS))))以上のコードは以下のコードでも同じく作動します。
ipw_df <- ipw_df %>%
mutate(W = treat * (1 / PS) + (1 - treat) * (1 / (1 - PS)))むろん、dplyrのmutate()を使わずに、以下のように入力しても重み変数は作成可能である。
ipw_df$W <- ifelse(ipw_df$treat == 1,
1 / ipw_df$PS,
1 / (1 - ipw_df$PS))最後に処置効果の推定です。CEMの時と同じやり方です。weight = ...引数を追加するだけです。
# 回帰分析にweight = ...を指定するだけ
IPW_fit <- lm(re78 ~ treat, data = ipw_df, weights = W)
modelsummary(IPW_fit)| Model 1 | |
|---|---|
| (Intercept) | 6422.839 |
| (397.427) | |
| treat | 224.676 |
| (577.657) | |
| Num.Obs. | 614 |
| R2 | 0.000 |
| R2 Adj. | -0.001 |
| AIC | 12796.0 |
| BIC | 12809.3 |
| Log.Lik. | -6394.998 |
| F | 0.151 |
推定された処置効果は約224.676ドルです。
それでは{WeightIt}パッケージを使ってみましょう。これまで使ってきました{MatchIt}パッケージや傾向スコア算出のためのglm()関数と使い方が非常に似ています。まず、第一引数として処置変数を結果変数、処置有無に影響を与えると考えられる共変量を説明変数とした数式オブジェクト(formulaクラス)を指定します。続いて、データ(data)、IPW算出の方法(method)、推定の対象(estimand)を指定します。データはdata = lalonde_dfとし、傾向スコアからIPWを算出するためmethod = "ps"を指定、最後にATT推定のためにestimand = "ATT"を指定します。ATEあるいはATCを計算する場合は"ATT"を適宜変更してください。また、IPW算出のために今回は傾向スコアを使いましたが、「処置を受ける確率」が計算できるなら何でも良いです。たとえば、Imai and Ratkovic (2014)10が推奨しているCovariate Balancing Propensity Score (CBPS) を使用する場合は"ps"の代わりに"cbps"を、複数の推定を組み合わせるスーパーラーニングをする場合は"super"11などが使えます。他にもエントロピーバランシングなど様々なオプションが提供されています。ここでは、先程の手法との比較のために、デフォルトの"ps"を使います。
Weighted_Data <- weightit(treat ~ age + educ + black + hispanic + married +
nodegree + re74 + re75, data = lalonde_df,
method = "ps", estimand = "ATT"){WeightIt}パッケージの便利なところはcobalt::love.plot()によるバランスチェックが簡単ということです。実際にやってみましょう。
love.plot(Weighted_Data, thresholds = 0.25, abs = TRUE)最後にIPW推定量を計算してみましょう。ここは先ほどの面倒くさい方法と同じですが、重み付け変数はweightit()関数から得られたオブジェクト(Weighted_Data)のweights列を使用します。推定ができたら結果を出力します。
IPW_Result <- lm(re78 ~ treat, data = lalonde_df,
weights = Weighted_Data$weights)
modelsummary(IPW_Result)| Model 1 | |
|---|---|
| (Intercept) | 5135.072 |
| (401.194) | |
| treat | 1214.071 |
| (568.904) | |
| Num.Obs. | 614 |
| R2 | 0.007 |
| R2 Adj. | 0.006 |
| AIC | 13241.9 |
| BIC | 13255.1 |
| Log.Lik. | -6617.942 |
| F | 4.554 |
という結果が得られましたが、これは{WeightIt}パッケージを使わなかった結果と完全に一致しますね。
ここではマハラノビス最近傍マッチングの例を使って、ATTを推定してみましょう。まずは、処置群の潜在的結果を統制群から割り当てて差分を算出し、その平均値を求るというのがノンパラメトリックな推定方法です。以下のコードが理解できるようになれば、他の最近傍マッチング(傾向スコアマッチングなど)も外部パッケージを使わずに計算できるようになります。
各処置群がどのケースとマッチングされたかはマッチングのオブジェクト名の後に$match.matrixを付けることで抽出できます。ちょっとやってみましょう。
# 全部表示させると長くなんるので、最初の6行まで確認
head(MH_Matching$match.matrix) ## [,1]
## 1 "553"
## 2 "557"
## 3 "478"
## 4 "567"
## 5 "579"
## 6 "438"
左が処置群、右が統制群ですね。それぞれケースのIDを抽出しましょう。
# MH.Matching$match.matrixは1列の行列です。
# 処置群IDは行の名前として格納され、統制群IDは1列目に入っています。
MH_Treat <- as.vector(row.names(MH_Matching$match.matrix))
MH_Control <- as.vector(MH_Matching$match.matrix[, 1])
# 中身を確認
print(MH_Treat)## [1] "1" "2" "3" "4" "5" "6" "7" "8" "9" "10" "11" "12"
## [13] "13" "14" "15" "16" "17" "18" "19" "20" "21" "22" "23" "24"
## [25] "25" "26" "27" "28" "29" "30" "31" "32" "33" "34" "35" "36"
## [37] "37" "38" "39" "40" "41" "42" "43" "44" "45" "46" "47" "48"
## [49] "49" "50" "51" "52" "53" "54" "55" "56" "57" "58" "59" "60"
## [61] "61" "62" "63" "64" "65" "66" "67" "68" "69" "70" "71" "72"
## [73] "73" "74" "75" "76" "77" "78" "79" "80" "81" "82" "83" "84"
## [85] "85" "86" "87" "88" "89" "90" "91" "92" "93" "94" "95" "96"
## [97] "97" "98" "99" "100" "101" "102" "103" "104" "105" "106" "107" "108"
## [109] "109" "110" "111" "112" "113" "114" "115" "116" "117" "118" "119" "120"
## [121] "121" "122" "123" "124" "125" "126" "127" "128" "129" "130" "131" "132"
## [133] "133" "134" "135" "136" "137" "138" "139" "140" "141" "142" "143" "144"
## [145] "145" "146" "147" "148" "149" "150" "151" "152" "153" "154" "155" "156"
## [157] "157" "158" "159" "160" "161" "162" "163" "164" "165" "166" "167" "168"
## [169] "169" "170" "171" "172" "173" "174" "175" "176" "177" "178" "179" "180"
## [181] "181" "182" "183" "184" "185"
print(MH_Control)## [1] "553" "557" "478" "567" "579" "438" "592" "528" "590" "564" "540" "566"
## [13] "561" "608" "463" "552" "556" "544" "578" "547" "460" "473" "585" "537"
## [25] "596" "563" "603" "524" "577" "525" "586" "604" "516" "558" "454" "581"
## [37] "523" "571" "457" "527" "568" "422" "583" "450" "573" "359" "449" "451"
## [49] "569" "570" "532" "565" "559" "374" "601" "594" "411" "600" "575" "538"
## [61] "550" "420" "382" "384" "613" "468" "421" "464" "539" "584" "477" "588"
## [73] "599" "430" "415" "396" "533" "526" "409" "574" "377" "426" "546" "609"
## [85] "447" "529" "397" "508" "406" "509" "506" "598" "610" "522" "560" "513"
## [97] "555" "593" "365" "429" "462" "458" "459" "542" "424" "511" "452" "562"
## [109] "400" "474" "504" "576" "520" "456" "445" "507" "344" "391" "443" "366"
## [121] "405" "461" "387" "497" "475" "351" "501" "370" "363" "399" "436" "311"
## [133] "362" "317" "352" "340" "410" "390" "353" "280" "442" "302" "333" "304"
## [145] "329" "321" "325" "338" "303" "319" "275" "301" "327" "328" "323" "339"
## [157] "343" "314" "259" "265" "320" "324" "277" "291" "257" "232" "271" "255"
## [169] "251" "247" "274" "261" "217" "235" "243" "220" "233" "269" "216" "202"
## [181] "192" "208" "191" "186" "193"
次は、元のlalondeデータから処置群と統制群の結果変数 (re78)を抽出し、MH.Matched.dfというデータフレームとしてまとめます。続いて、各ケースのITEを計算します。
MH_Matched_df <- tibble(
Treat_ID = MH_Treat,
Control_ID = MH_Control,
Treat_Y = lalonde[MH_Treat,]$re78,
Control_Y = lalonde[MH_Control,]$re78
)
MH_Matched_df <- MH_Matched_df %>%
mutate(ITE = Treat_Y - Control_Y)MH_Matched_df1行目を見ると1は553とマッチングされ、それぞれの結果変数は9930.046ドル、0ドルです。その差分がITEであり、9930.046ドルですね。
これらITEの平均値がATTとなります。早速計算してみましょう。
MH_ATT <- mean(MH_Matched_df$ITE)
MH_ATT## [1] 507.0434
ATTは約507.043ドルでした。回帰分析でやった結果と一致しますね。大人しく回帰分析でやりましょう。
これまでは主にQ-Q Plotを用いたバランスチェックを紹介しましたが、MatchItパッケージは標準化差分/標準化バイアスのバランスチェックにも対応します。やり方としては、
X <- matchit(...)
X_Blc <- summary(X, standardize = TRUE)
plot(X, interactive = FALSE)です。たとえば、傾向スコアマッチングの例だと、
PS_Blc <- summary(PS_Matching, standardize = TRUE)
plot(PS_Blc, interactive = FALSE)
abline(v = 0.25, col = "red")interactive = TRUEにすると各ポイントがどの変数に対応するかが分かります。灰色のポイントはバランスが改善された共変量であり、黒は悪化した共変量です。また、これは前の講義で解説した標準化差分ですが、100を掛けていないものであり、なお、絶対値です。したがって、バランスの基準は0.1となります。多くの共変量においてバランスが改善が見られますが、0.1未満の共変量も少なく、大きな改善とは考えられませんね。
もっと分かりやすい図がほしければ、これまで何回か使ってきましたcobaltパッケージがおすすめです。
love.plot(PS_Matching, threshold = 0.25, abs = TRUE)赤い点がマッチング前の標準化差分、青い点がマッチング後の標準化差分です。図の上にあるdistanceは傾向スコアを意味します。educとhispanはマッチングによってバランスが悪化しましたが、他の変数が改善された幅を考えると許容範囲かもしれません。それでも、0.25の垂直線の右側にある青い点が3つもあり、傾向スコアの標準化差分もかなり大きいのが分かります。
他にも傾向スコアをマッチング前後に分けてヒストグラム化することもいいでしょう。
bal.plot(PS_Matching, var.name = "distance", which = "both",
type = "histogram", mirror = TRUE)左がマッチング前、右が後です。また、赤いヒストグラムは統制群、青は処置群です。もし、傾向スコアのバランスが取れているならヒストグラムは上下に対称となります。マッチング後の傾向スコアの標準化差分はかなり大きかったのですが、それでも大幅の改善は見られました。少なくとも、普通に回帰分析するよりも良さそうですね。
これまでの推定値を同時に示す図を作ってみましょう。そのためには各推定結果のオブジェクトから処置変数の係数を抽出する必要があります。たとえば、重回帰分析によるATEの推定結果はRaw_Fit2に格納されています。lm()やglm()関数で行った推定結果から係数のみを抽出する関数はcoef()です。
coef(Raw_Fit2)## (Intercept) treat age educ black
## 66.5145197 1548.2438020 12.9776337 403.9412309 -1240.6440820
## hispanic married nodegree re74 re75
## 498.8968650 406.6208454 259.8173680 0.2963774 0.2315259
各係数の推定値がベクトルとして表示されます。ここで処置変数(treat)の係数はベクトルの2番目の要素であるため、coef(Raw_Fit2)[2]で抽出できます。他の手法による推定についても同じです。なぜなら、これまで処置変数を1番目の説明変数として使ってきたからです。
一方、標準誤差の抽出はちょっと面倒くさいです。いくつかの方法がありますが、ここではsummary(オブジェクト名)$coefを使います。推定結果を出力するsummary()関数のcoef要素を抽出すると以下のようなデータフレームが抽出されます。
summary(Raw_Fit2)$coef| Estimate | Std. Error | t value | Pr(>|t|) | |
|---|---|---|---|---|
| (Intercept) | 66.5145197 | 2436.7458028 | 0.0272965 | 0.9782323 |
| treat | 1548.2438020 | 781.2792975 | 1.9816778 | 0.0479682 |
| age | 12.9776337 | 32.4889061 | 0.3994482 | 0.6897042 |
| educ | 403.9412309 | 158.9062381 | 2.5420099 | 0.0112704 |
| black | -1240.6440820 | 768.7644071 | -1.6138157 | 0.1070897 |
| hispanic | 498.8968650 | 941.9425489 | 0.5296468 | 0.5965515 |
| married | 406.6208454 | 695.4723228 | 0.5846686 | 0.5589889 |
| nodegree | 259.8173680 | 847.4420524 | 0.3065901 | 0.7592610 |
| re74 | 0.2963774 | 0.0582726 | 5.0860483 | 0.0000005 |
| re75 | 0.2315259 | 0.1046199 | 2.2130198 | 0.0272694 |
ここで、treatの標準誤差は2行2列目に格納されていることが分かります。したがって、treat標準誤差はsummary(Raw_Fit2)$coef[2, 2]で抽出できます。他にもsummary(Raw_Fit2)から抽出できる要素は色々あります。興味のある方はstr(summary(Raw_Fit2))で確認してみてください。
Comparison_df <- data.frame(
Method = c("Single Regression",
"Multiple Regression",
"Mahalanobis (ATT)",
"CEM (ATT)",
"Propensity Score (ATT)",
"IPW (ATT)"),
Est = c(coef(Raw_Fit1)[2],
coef(Raw_Fit2)[2],
coef(MH_Fit)[2],
coef(CEM_Fit)[2],
coef(PS_Fit)[2],
coef(IPW_fit)[2]),
SE = c(summary(Raw_Fit1)$coef[2, 2],
summary(Raw_Fit2)$coef[2, 2],
summary(MH_Fit)$coef[2, 2],
summary(CEM_Fit)$coef[2, 2],
summary(PS_Fit)$coef[2, 2],
summary(IPW_fit)$coef[2, 2])
)続いて、各手法名(Method列)を順序付きの変数に変換します。これをしないと、各手法がアルファベット順で並ぶことになります。順番が、データフレームに登場する順番の逆となります。したがって、fct_inorder()を使い、「登場順でfactor化」をし、fct_rev()で順番の反転を行います。fct_*()関数群はforcatsパッケージが提供しており、詳細は宋・矢内の第14章を参照してください。
Comparison_df <- Comparison_df %>%
mutate(Method = fct_inorder(Method),
Method = fct_rev(Method))続きまして、95%信頼区間の下限(LL)と上限(UL)を追加し、\(\alpha = 0.05\)水準で統計的有意か否かを示す変数も作成します。
Comparison_df <- Comparison_df %>%
mutate(LL = Est + qnorm(0.025) * SE,
UL = Est + qnorm(0.975) * SE,
Sig = if_else(LL * UL > 0 , "Significant", "Insignificant"))Comparison_df| Method | Est | SE | LL | UL | Sig |
|---|---|---|---|---|---|
| Single Regression | -635.0262 | 657.1374 | -1922.99190 | 652.9395 | Insignificant |
| Multiple Regression | 1548.2438 | 781.2793 | 16.96452 | 3079.5231 | Significant |
| Mahalanobis (ATT) | 507.0434 | 746.8731 | -956.80109 | 1970.8878 | Insignificant |
| CEM (ATT) | 1070.9073 | 1248.1293 | -1375.38120 | 3517.1959 | Insignificant |
| Propensity Score (ATT) | 894.3675 | 730.2948 | -536.98398 | 2325.7189 | Insignificant |
| IPW (ATT) | 224.6763 | 577.6575 | -907.51150 | 1356.8641 | Insignificant |
宋のように「表より図がいい!」という方は以下のように可視化することも可能です。
Comparison_df %>%
ggplot() +
# x = 0の垂直線を引く。線の種類は破線 (linetype = 2)で
geom_vline(xintercept = 0, linetype = 2) +
# 点と線を同時に出力するgeom_pointrange幾何オブジェクトを使用
# 点の位置 (x, y)はEst, Methodに、線の下限と上限 (xmin, xmax)はLL, UL
# Sigの値によって色分けするためにcolor引数もマッピング
geom_pointrange(aes(x = Est, y = Method, xmin = LL, xmax = UL,
color = Sig), size = 0.75) +
# colorにマッピングされているSigの値がSignificantなら色をblackに
# Insignificantならgray70にする
scale_color_manual(values = c("Significant" = "black",
"Insignificant" = "gray70")) +
# ラベルを付ける
labs(y = "Methods", x = "Estimates (ATE or ATT)",
color = "") +
# minimalテーマを適用
theme_minimal() +
# 文字の大きさを調整
theme(text = element_text(size = 16))lalondeデータセットを提供するパッケージは複数あり、それぞれデータセットの構成に違いがあります。どれを使っても実習には問題ございませんが、以下の内容を再現される場合はdata()内にpackage = "cobalt"を指定するようにしましょう。↩︎
ここでは基準としている標準化差分は0.25ですが、{BalanceR}の25と同じです。↩︎
これはmatchit()内にreplacement = TRUEを指定することで防ぐことができます。既定値はFALSEですが、TRUEを指定すれば、マッチングされた統制群ケースを2回以上マッチングすることができます。↩︎
これはmatchit()内にratio引数を指定することで防ぐことができます。既定値は1であり、この場合は1:1マッチングです。↩︎
これは重み変数を使わない最近傍マッチングをしても追加されます。ただし、その場合、weightsは全て1になります。↩︎
実はnearest = "logit"が省略されています。つまり、ロジスティック回帰分析から得られた傾向スコアの距離に基づくマッチングを意味します。↩︎
Austin, Peter C., and Guy Cafri. 2020. “Variance Estimation When Using Propensity-Score Matching with Replacement with Survival or Time-to-Event Outcomes.” Statistics in Medicine, 39 (11): 1623–40.↩︎
ただし、replacement = TRUEを指定した場合、重み付け回帰分析をする必要があります。。なぜなら、統制群から1つのケースが何回もマッチングされる可能性があるからです。統制群のあるケースが数回マッチングされてもPS_Dataにそのケースは1つしかありません。この場合は重みをより高く付ける必要があります。幸い、{MatchIt}パッケージはその重みの自動的に計算してくれます。具体的にはCEMと同様、weights変数がPS_Dataに追加されます。↩︎
ロジスティック回帰分析ではなくても問題ありません。他、プロッビット、サポートベクトルマシン、回帰木、判別分析、ニューラル・ネットワークなどもあり、これらを組わせることも可能です。↩︎
Imai, Kosuke and Marc Ratkovic. 2014. “Covariate Balancing Propensity Score.” Journal of the Royal Statistical Society, Series B, Vol. 76, No. 1, pp. 243-246.↩︎
Super Learnerを使った例は「Cyrus, Samii, Laura Paler, and Sarah Zukerman Daly. 2016. “Retrospective Causal Inference with Machine Learning Ensembles: An Application to Anti-recidivism Policies in Colombia.” Political Analysis, 22 (4) pp. 434-456」を、日本語による解説は誰かの報告スライドを参照して下さい。↩︎